// deduce flux from equation and finite-volume alogrithm
namespace Dune {
  namespace Flux {
    template<typename EQ, typename NUMBER>
    class FluxDeduce {
    public:
      using RangeField = NUMBER;
      using RF = RangeField;
      using DF = RF;
      static constexpr int dim = EQ::dim;
      static constexpr int nvar= EQ::nvar;
      using Range = Dune::FieldVector<RangeField,nvar>;
      FluxDeduce(EQ& eq, Dune::ParameterTree& ptree) {
        this->eq = &eq;
        this->ptree = ptree;
      }
      virtual ~FluxDeduce() {
      }
  
      template<typename E, typename X>
      void algoFlux(const E& inside, const X& x_inside,
                    const E& outside, const X& x_outside,
                    const Dune::FieldVector<DF,dim> n_F,
                    const Dune::FieldVector<RF,nvar>& u_s,
                    const Dune::FieldVector<RF,nvar>& u_n,
                    Dune::FieldVector<RF,nvar>& f) const
      {
        // fetch flux
        Dune::FieldMatrix<RF,nvar,dim> Fs;
        Dune::FieldMatrix<RF,nvar,dim> Fn;
        // evaluate flux
        equFlux(inside,x_inside,u_s,Fs);
        equFlux(outside,x_outside,u_n,Fn);

        // max eigenvalue
        RF alpha_s(0.0);
        RF alpha_n(0.0);
        double u_s_tmp[nvar];
        double u_n_tmp[nvar];
        for (int i=0; i<nvar; i++) {
          u_s_tmp[i] = u_s[i];
          u_n_tmp[i] = u_n[i];
        }
        alpha_s = eq->local_max_eigenvalue(u_s_tmp);
        alpha_n = eq->local_max_eigenvalue(u_n_tmp);

        f = 0.0;
        // scheme: lax-fredrich, local lax-fredrich
        std::string algofluxname = ptree.get<std::string>("fem.algofluxname");
        double cfl = ptree.get<double>("fem.cfl");
        double dx = ptree.get<double>("grid.structured.LX") 
          / ptree.get<int>("grid.structured.NX"); 
        if (algofluxname == "lf" || algofluxname == "llf") { 
          Fs.umv(n_F,f);
          Fn.umv(n_F,f);
          f *= 0.5;
          if (algofluxname == "lf")
            for (size_t i = 0; i<nvar; i++) {
              // f[i] = f[i] + dx / (dim * dt) / (2.0 * dim) * (u_s[i] - u_n[i]);
              f[i] = f[i] + 2.0 * 0.5 * (u_s[i] - u_n[i]);
            }
          if (algofluxname == "llf")
            for (size_t i = 0; i<nvar; i++)
              f[i] = f[i] + std::max(alpha_s,alpha_n) / cfl / (2.0 * dim) * (u_s[i] - u_n[i]);
        }
      }
 
      template<typename E, typename X, typename RF>
      void equFlux (const E& e, const X& x,
               const Dune::FieldVector<RF,nvar>& u,
               Dune::FieldMatrix<RF,nvar,dim>& F) const
      {
        double utmp[nvar],ftmp[nvar];
        for (int i=0; i<nvar; i++) {
          utmp[i] = u[i];
          ftmp[i] = 0.0;
        }
        for (int j=0; j<dim; j++) {
          eq->get_equation_flux(j, ftmp, utmp);
          for (int i=0; i<nvar; i++)
            F[i][j] = ftmp[i];
        }
      }

      void set_dt(RF dt) {
        this->dt = dt;
      }

    private:
      Dune::Equation::EquationBase* eq;
      Dune::ParameterTree ptree;
      double dt;
    };
  }
}
